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Abstract 


This  paper  discusses  a procedure  for  interval  estimation  of  the 
mean  0 of  a correlated  binary  (0,1)  sequence.  The  method  assumes  that 
the  sequence  is  strictly  stationary  and  that  a particular  string  of  m 
binary  digits  is  a recurrent  event  in  the  sequence,  where  nn>l  is  unknown. 

Of  the  2m  choices  for  the  possible  recurrent  events,  the  strings  of  all 
zeros  and  of  all  ones  are  examined. 

For  each  m=l,2,...  the  sequence  is  demarcated  by  entrance  to  the 
recurrent  event.  The  subsequences  between  the  demarcation  points  thus 
form  independent  epochs  by  assumption.  Classical  techniques  then  yield, 
variance  estimates  for  the  number  of  ones  and  zeros  in  the  epoch  as  well 
as  an  estimate  of  the  covariance  of  the  ones  and  zeros.  A quadratic 
equation  in  0 is  solved  to  obtain  an  interval  estimate. 

Each  string  of  all  ones  or  all  zeros  examined  yields  a 1-a  confi- 
dence interval.  The  intervals  are  intersected  to  obtain  shorter  intervals 
with  confidence  greater  than  l-2a.  Since  each  m=l,2,...  yields  an  interval, 
a conservative  rule  is  developed  to  determine  the  m whose  interval  is 
finally  used.  This  rule  is  based  upon  the  empirical  run  lengths  in  the 
binary  sequence. 

The  procedure  is  then  applied  to  interval  estimation  of  the  fractile 
for  the  waiting  time  distribution  in  a simulation  of  the  M/M/1  queue  with 
activity  level  0.9.  For  0=0.1  and  0.5  the  proposed  method  worked  well. 

For  0=0.9  results  showed  some  degradation.  An  error  analysis  led  to  a 
set  of  recv Timendations  for  keeping  performance  in  practice  close  to  the 
desired  theoreti',i  level:.  An  appendix  describes  algorithms  for  computing 
the  critical  quantities  upon  which  the  proposed  method  relies. 
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1 . Introduction 

Although  many  techniques  exist  for  the  statistical  analysis  of 
simulation  output  in  general  [11],  little  has  appeared  in  the  literature 
that  specifically  addresses  the  estimation  of  fractiles  of  a distribution 
that  arises  in  a simulation.  This  problem  is  not  unique  to  simulation. 

In  its  more  general  form  it  concerns  the  estimation  of  the  mean  of  a 
random  binary  (0,1)  sequence  whose  elements  are  correlated.  The  purpose 
of  this  paper  is  to  describe  a method  of  interval  estimation  for  the  more 
general  problem  and  illustrate  the  method  with  a simulation  example. 

The  proposed  method  relies  on  the  theory  of  recurrent  events,  as 
described  in  Feller  [9].  Although  application  of  these  methods  to 
simulation  output  analysis  is  not  new,  the  approach  taken  here  is. 

Crane  and  Iqlehart  [7,8]  and  Fishman  [11]  exploit  the  regenerative  prop- 
erties of  certain  simulation  models  to  analyze  output.  The  theory  of 
regenerative  processes  generalizes  the  theory  of  recurrent  events  to  in- 
clude continuous  as  well  as  discrete  time.  See  Smith  [13]  for  details. 

What  is  new  here  is  the  application  of  the  theory  of  recurrent  events 
without  reference  to  particular  simulation  models.  Specifically,  the 
proposed  method  enables  one  empirically  to  identify  a state  that  appears 
to  have  the  nrrmorties  of  a recurrent  state  and  then  to  use  that  state  to 
cut  up  the  sample  path  into  approximately  independent  identically  distri- 
buted segments  that  can  then  be  used  with  relatively  elementary  statistical 
methods  to  compute  an  interval  estimate.  Moreover,  the  binary  nature  of 
the  data  offers  considerable  computational  conveniences  that  make  this 
method  of  inference  attractive. 

Whenever  one  deals  with  discrete  finite  state  data,  the  inclination 
to  approximate  the  sequence  by  a Harkov  chain  is  hard  to  resist.  Once  the 
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order  of  the  chain  is  determined,  the  subsequent  inference  is  well  known  [1]. 
Moreover,  every  state  of  the  given  order  is  a recurrent  state.  To  determine 
the  order  of  the  chain  one  can  apply  chi-square  or  likelihood  ratio  methods 
the  theory  for  which  is  also  known  in  principle  [1].  Although  our  study 
began  in  this  way,  we  quickly  learned  that  the  computational  problems  that 
arose  in  estimating  the  order  were  considerably  greater  in  complexity  than 
those  of  the  less  restrictive  recurrent  event  approach.  Moreover,  the 
performance  was  not  nearly  as  good  as  with  the  recurrent  event  approach. 
Therefore  we  present  here  the  results  for  the  recurrent  approach,  even 
though  the  statistical  inference  for  this  theory  is  not  as  complete  in  the 
statistical  literature  as  for  Markov  chains. 

Section  2 introduces  the  reader  to  the  problem  as  it  is  formulated  in 
the  theory  of  recurrent  events.  Section  3 describes  how  one  can  obtain 
shorter  interval  estimates  by  using  results  based  on  analyses  for  different 
recurrent  states.  Section  4 indicates  which  among  the  many  potential 
states  deserve  consideration.  Section  5 describes  estimators  that  need  to 
be  substituted  for  population  parameters  and  Section  6 contains  a computing 
scheme  that  facilitates  efficient  computation.  Section  7 presents  r^s"lts 
that  show  how  the  method  performs  for  the  .1,  .5,  and  .9  fractiles  of  the 
waiting  time  distribution  in  an  M/M/1  queueing  simulation  with  activity 
level  of  .9.  Finally,  Section  8 contains  recommendations  for  inter- 
pretation and  utilization  of  interval  estimates  that  a user  of  the  pro- 
posed method  may  encounter  in  practice. 

?.  The  Recurrent  Model 

Let  {X.;  I , . . . ,n , denote  a sequence  of  observations  on  a 
strictly  stationary  binary  (0,1)  stochastic  sequence.  Define  for  m n 
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v!m)  = mf1  2jX 

1 j=0  1-J 


i -m j » • • f n • 


Suppose  there  exists  an  m and  0<y  <2-1  such  that 

(2)  pr(Xk  = x | Xj  = z,  Y|m)=  y*)  = pr(Xk  = x I y*) 


j<i<k;  n>i>m; 


x,z  = 0,1 . 


One  can  then  apply  the  theory  of  recurrent  processes  [13]  to  the  analysis 

(m  J * 

of  {X^},  provided  Y:  = y infinitely  often  as  n -*■  °°.  In  particular. 


«jm)  h j 5(v!m)  - y*) 
1 J 


i=m, . . . ,n 


6(X)  5 


1 x=0 

0 x*0 


lSm)  5 min  film:  = j) 

J 


j = l , . . . ,n' 


1,  T(7L  = n+l. 

%W 


Here  {T\m'  ; j^l , . . . ,N'm^ } defines  a discrete  renewal  process  where 
J n 

is  the  time  of  the  jth  renewal. 

J 


c(m)  _ T(m)  T(m) 

j J+T  J 


Am)  _ f x 

J i^fm)  Xi 


Tfm)  , 

j+l 


j=0,...,N 


Then  {C\m)  ; j=l ....  ,N^m,_i } and  fS  ; j=l N^-i}  are  each 

sequences  of  i.i.d.  non  negative  integer  valued  random  variables.  Let 
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(11) 


N -1 
n 


N -1 
n 1 


C s l C. 
j=l  J 


l Sj 

j = l 3 * 


Q being  the  l-a/2  value  of  the  unit  normal  distribution.  Then  for  large 
n,  pr[f(ohO]  4 1 -a . Solutions  to  the  probability  argument  yield  the  seven 
cases  in  Table  1.  Here  a1  and  a^  are  the  roots  of  f(e) 


lable  1 

I -a  Interval  Estimate  of  0 


Case 

Interval 

f (0) 

f(0) 

fll) 

f'Ol- 

1 

[a-,  ,a2] 

+ 

- 

+ 

+ 

2 

[a i ,1] 

+ 

- 

- 

+ 

+ 

3 

[0,a2] 

- 

+ 

+ 

4 

[0,a1]u[a2,l] 

- 

+ 

- 

- 

+ 

5 

[0,a^] 

- 

+ 

+ 

6 

[a2,l] 

+ 

+ 

- 

- 

7 

[0,1] 

f (©)  f o vo<  [o.i] 

with 


'1 


<_  a2  • One  can  easily  show  that  as  n ->  ® case  1 becomes  the 


type  of  interval  estimate  obtained.  This  procedure  for  computing  interval 
estimates  is  due  to  an  extension  of  the  work  in  Bliss  [2,3]  by  Fieller[10]. 


i 


3.  Many  Recurrent  States 

In  practice  many  stochastic  processes  with  a discrete  state  space 
have  more  than  one  recurrent  state.  In  particular,  if  the  sequence  of 
states  {y j = j-1;  j=l,...,2m}  were  all  recurrent  then  {X^l  would  be 

an  mth  order  Markov  chain.  Although  an  analyst  is  free  to  choose  any 
of  the  set  of  recurrent  states  for  the  computation  of  an  interval  estimate 
each  state  produces  an  interval  of  different  width.  Moreover,  the  n for 
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which  the  aforementioned  asymptotic  properties  usually  hold  differs  with 
each  selected  state. 

Suppose  that  for  a given  {X.}  one  computes  a set  of  k interval 
estimates  for  o.  For  expository  convenience  assume  that  all  the  intervals 
correspond  to  case  1 in  Table  1 and  that  [a^ (i ) ^(i )]  is  the  ith 
1-a  interval  estimate  of  0 for  i=l,...,k.  Then,  using  Bonferroni 's 
inequality  [12]  one  can  show  that 

k 

(12)  n pr[a,(i)  < 0<  ap( i )]  ^ 1-ka  . 

i = l 1 

If  these  intervals  intersect  then  the  intersection  provides  an  interval 
estimate  for  o with  probability  exceeding  1-ka. 

Consider  the  case  k=2  for  which  a^ (1 )sa^ \ 2 ) ^a^ ( • )sa^(2) . Then 
[a, (2) ,a2(l )]  includes  e with  probability  of  at  least  l-2a  . Here 
a question  arises  as  to  whether  a shorter  interval  would  be  obtained  by  computing 
a l-2a  interval  estimate  for  i=l  or  i=2  separately.  Section  8 ex- 
amines this  issue  with  regard  to  a specific  example. 

4.  States  to  Consider 

As  the  introduction  indicates,  our  intention  is  to  describe  a procedure 
for  computing  interval  estimates  for  0 by  approximating  ( Xj > by  a re- 
current process.  Since  the  state  vector  as  defined  in  (1)  implies  2 potential 
recurrent  states,  even  a moderate  m yields  an  excessive  number  of  options. 

it 

However,  two  particular  states  deserve  special  attention.  I he  state  y = 0 
implies  t it  when  Y.=y!  X.  .=0  for  j=l,...,m.  Now  it  is  plausible  that 
for  sufficient!'  .arge  p Y^y*  contains  all  conditioning  information 
about  the  past  and  therefore  this  y*  of  order  m enables  one  to 
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demarcate  { X ^ } into  sequences  of  i.i.d  random  variables  as  in  (4). 

"ft  m ★ 

Alternatively  the  state  y = 2-1  implies  that  when  Y^y,  ^ = 1 for 

j=l , . . . ,m.  A similar  plausibility  argument  applies  here. 

in  the  remainder  of  this  paper  we  restrict  attention  to  these  poten- 
tial regenerative  states.  Doing  so  relieves  us  of  the  burden  of  consid- 
ering the  remaining  2m-3  states  for  each  value  of  m considered  but 
leaves  us  open  to  the  possibility  of  excessive  computation  if  m turns 
out  to  be  large.  Happily  the  binary  character  of  the  data  makes  this 
possibility  remote  as  the  computing  scheme  in  Section  6 srows. 


5.  Sample  Variances 

In  practice  one  does  not  know  t(N)»  0 , o,^  and  o^.  To  resolve  this 

problem  we  consider  the  sample  quantities 

N -1 

s2(c)  = i (c.  -c r 
j=i  j 

N -1 

o h p 

S2(s)  = l (S.  -S)d 

j=l  J 

Nn-1 

s(C,S)  = l (c • - C) (Si  - s) 

0=1  J J 

N -1 

' ' fa  ^ ’ 


(13) 


Since 


S = 


V1 


V1 


N -1 
n 


N -1 
n 


N -1 
n 


One  can  use  the  well  known  results  of  renewal  theory  [5,13] 


E<V  ~ " / 


(14a) 
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(14b) 


var(Nn) 


"°CC 


in  a series  expansion  to  show  that 

(15)  E[s2(C)]  = E(Nn-l)occ  - 2occ  + Oll/n). 

Similarly 

E[s2(S)]  = E(Nn-l)oss  - as$  - e2acc  + 0(l/n) 

(16) 

E[s(C,S)3  = t(Nn-l)oc$  - ocs  - 0acc  + 0(l/n). 

These  results  reveal  that  s2(C),  s2(S)  and  s(C,S)  are  asymptotical ly 

unbiased  estimations  of  E(Nn-l)occ>  E(Nn_1^°SS  and  ^Nn_1)0CS  re~ 
spectively.  Moreover,  one  can  substitute  these  quantities  Into  (9b)  and 
show  that  the  resulting  distribution  of  Z(m)//var  Z(ml  again  approaches 

the  standard  normal  law  as  m °°  . 


t 

6.  A Computing  Scheme 

For  analysis,  binary  data  offer  conveniences  that  allow  considerable 
computational  efficiency.  The  sequence  {X..}  consists  of  runs  of 
ones  and  zeros  of  varying  lengths 

i length  of  the  jth  run  of  zeros 

(17)  J' 

L2j  = Length  of  the  jth  run  of  ones 

where  L]  = 0 if  X^l.  If  M runs  occur  then  {L1 ,LM)  summarizes 
{X.)  without  any  information  loss.  One  can  easily  establish  a bound  on 
E(M).  Suppos  {X.}  i a sequence  of  i.i.d  random  variables  with 
prlx.=l)=p.  Then  E(M)  = np(l-p)  < n/4.  If  p=.9  E(M)  = 0.09n. 

+ The  Appendix  contains  algorithms  for  computing  the  critical  quantities 
described  here. 
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Now  in  the  cases  at  hand  we  anticipate  positive  correlation  between 
elements  of  {X^.  Therefore  one  expects  runs  to  be  longer  and 
be  smaller  for  a given  n. 

The  computing  scheme  presented  here  uses  the  recurrent  states 

"k  J/m  , 

y = 2 -1  for  k = 0 and  1.  In  the  remainder  of  this  section  we  assume 
k and  m given  and  suppress  reference  to  them  except  where  clarity 
calls  for  explicit  mention.  Let 


. k - 1 


Kj  = min  (2i  + k-1  > K^:  ^-j+k-i  - m)  J = 1»2»' 

For  computational  convenience  we  set  k.  s M+l  if  the  aforementioned 

J 


minimum 


does  not  exist.  Define 


(19)  J = min  (j  : Kj  = M+l  ) - 1 

Then  (K.;  j=l,...,J}  forms  a sequence  of  cutpoints  for  given  m and 
k.  As  m increases  for  given  k the  number  of  cutpoints  decreases, 
which  accelerates  subsequent  search  procedures. 

For  given  k and  m let 


B,  = I Lk  - Jm 
J i=l  i 


Then  one  can  show 


J-l 


C = B 


(22) 


S = kB 


+ l D 

+ (l-k)‘  X Di  + (2k- 1 ) l E. 

J j=l  J 3=1  J 


J-l 


J-l 


N = B,  + J 
n J 


where  (11)  and  (3)  define  S,  C and  Nn-  One  can  also  show 

s2(C)  = B + Y D2  - C2/(N  -1) 

J j=l  J 


J-l 


(23) 


s‘(S)  = kBj  + J (kEj  + (1-k) (Dj-Ej ) } 
0 * 


- {kS2  + (1-k) (C-S)2}/(Nn-l ) 


J-l 

s(C,S)  = kBj  + l DjfkEj  + (1 -k) (Dj-Ej )} 
J ^ 


C{kS  + (l-k)(C-S)}/(Nn-l)  , 


where  (13)  defines  s2(C),  s2(S)  and  s(C,S). 


l,o  formulae  for  computing  these  quantities  apply  for  k = 0 and  1. 
Also  notic»  i-hat  the  number  of  calculations  J is  always  less  than  n 
and  usually  considerably  smaller.  An  additional  convenience  arises  for 

m < ^“2j+k-l  j“l  »•  • • iJ 


i (m , k ) ^ -j-^en  ^or  sc(,eme  ,n+i  and  k one  has 


j ( m+ 1 ,k)  _ j(m,k) 


,(m+l ,k)  _ „ f m , k ) 


,(m+l ,k)  . 


,(m+l  ,k)  _ 


(m+1 ,k) 


,(m+l  ,k)  _ 


_ j 

D(m,k)  + 1 
*3 

F(m,k)  , 
j 

N(m,k)  _ j(m,k) 


7.  Estimating  m 

The  foregoing  scheme  for  estimating  E(Nn-l)a££,  E(Nn-l)oj<.  and 
E(Nn-l)ocs  leaves  one  remaining,  but  critical,  problem  before  one  can 
compute  an  interval  estimate  for©.  This  concerns  the  estimation  of  m 
for  given  k.  Firstly,  m for  k=0  usually  differs  from  in  for  k=l . 
Secondly,  a criterion  is  needed  to  determine  a satisfactory  m for  each 
case.  Consider  the  linearized  form  (7)  with  variance  in  (9b).  Since 

(25)  s2(Z(m,k))  = s2(S(m,k))  - 20s(C(m,k,,S(m,k))  + §2s2(C'm,k)) 


orovides  an  estimate  of  (9b)  a conservative  rule  is 


(26)  m*  = min[m:  s2(Z(m’k)  2 s2(Z(i,k));  i=0,...,L*] 

Lk~:sup  lL2i  +k-l>  . 

ui^rM/21 

★ 

Here  denotes  the  longest  run  of  k's  in  the  sample  sequence  {L.l. 
This  rule  requires  total  enumeration  of  all  possible  run  lengths  to 


f The  quantity  fxl  denotes  the  integer  part  of  x. 
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determine  m*.  Moreover,  by  selecting  m*  based  on  the  maximal  s2(Z^i,k^) 
the  rule  picks  the  most  conservative  result  that  the  empirical  evidence 
will  support.  A less  costly  and  less  conservative  rule  is  obtained  as 
follows.  Let  m^  denote  the  ith  order  statistic  of  the  sequence 
^2i+k-l ,L2i+k-l  + i=l »• • • » C M/ 2 1 } . Then  the  rule  is  defined  as 

(27)  m**  = min  [m.:  s2(Z(mj’k))  ->  s2(Z(mi,k)),  i = l 2fM/2lj. 

J 

This  rule  compares  sample  variances  for  the  empirical  runs  of  k's  in 
the  sample.  These  are  a subset  of  the  runs  contained  in  (26). 

Regardless  of  which  rule  is  selected  to  estimate  m,  a possibility 
remains  that  the  sample  sequence  (L^>  does  not  contain  runs  of  sufficient 
length  to  estimate  m with  accuracy.  One  way  to  check  on  adequacy  is  to 
plot  s2(Z^,k^)  versus  i and  observe  if  tne  sample  variance  achieves 
an  approximate  plateau  in  the  vicinity  of  its  maximum.  If  it  does  then 
one  can  place  confidence  in  the  estimated  m as  providing  suitable  accuracy 
for  the  approximating  scheme. 

8.  A Queueing  Example 

To  illustrate  the  proposed  estimation  procedure  we  use  a simulation 
of  the  M/M/1  queueing  problem  with  arrival  rate  x = 1 and  service 
rate  u>  = 1/.9.  This  yields  an  activity  level  p = .9  which  implies  a high 
level  of  congestion.  In  particular  the  mean  number  of  jobs  in  queue  is 
p/(1-p)  = 9 and  the  mean  waiting  time  is  p/x ( 1 -p ) = 9.  The  objective  is  to 
estimate  fractiles  of  the  waiting  time  distribution  F corresponding  to 
waiting  times  w = 0,5.29,19.78.  These  fractiles  are  0 =.1,.5,.9,  respec- 
tively, which  cover  a substantial  part  of  the  domain  of  F. 

Let  (W.;  i = l,...,n}  denote  a sample  sequence  of  waiting  times  collected 
on  a given  simulation  run,  where  collection  beqins  after  the  efiects  of 

+ See  Cox  and  Smith  [6]. 
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initial  conditions  has  dissipated. 

(28)  y = I 1 


Defi ne 
W.  < w 

i 

W.  > w 

l 


so  that  0 in  (6)  estimates  the  fractile  of  F correspondino  to  w . 
Results  in  Blomquist  [4]  enable  one  to  show  that 


(29)  lim  n var(6)  = F*(w){l+p2-2F*(w)(l+A(l-p)w)/(l-p)2-|+F*(w) i 


n*» 


(30) 


F*lw)  = 


1 - P 

-wlw-A) 

pe 


w=0 

w>0 


so  that  one  can  compute  the  large  sample  variance  of  6 for  comparison 
with  empirical  results. 

For  each  value  of  0 the  sampling  experiment  consisted  of  n=8l92 
observations  on  100  independent  replication  collected  from  a simulation 
that  was  operating  in  the  steady  state.  Application  of  the  procedures 
described  in  Sections  5,6  and  7 produced  the  results  in  Table  2.  In 
particular,  the  theoretical  mean  is  known  a priori . The  sample  mean  denotes 
the  average  values  of  6 over  the  100  replications.  The  theoretical  variance 
follows  from  (29)  divided  by  n whereas  the  sample  variance  is  the 
average  value  of  s2(Z^n)  ’^)/(N^m  *^-l)n  over  the  100  repl  ications . 

The  coverage  rates  indicate  the  proportion  of  intervals  based  on  Table  1 
that  include  the  theoretical  o. 

Notice  that  the  coverage  rates  for  o = .1  and  .5  are  conservative 
whereas  those  for  o=.9  need  more  careful  scrutiny.  Although  the  sample 
variances  usually  underestimate  the  corresponding  theoretical  variance,  the 
width  of  the  intervals  appear  larger  than  theoretically  expected.  Ihe 
confusion  disappears  when  we  look  at  how  the  theoretical  interval  widths 
were  computed.  Asymptotically  one  can  show  that  on  a given  run  o has 


n 


95 

Order  of  : 
Schggie 


theoretical 

0.1000 

2.809 

0.95 

0.0657 

t 

n.a. 

Sample:  k=0 

0.1029 

1 .648 

0.99 

0.1269 

33 

k=l 

0.1029 

1 .937 

0.97 

0.1835 

l+t 

Theoretical 

0.5000 

34.00 

0.95 

0.2286 

n.a. 

Sample:  k=0; 

1 0.5016 

28.91 

0.98 

0.5050 

2b 

k=l 

0.5016 

28.48 

0.95 

0.4906 

2 

Theoretical 

0.9000 

29.54 

0.9b 

0.2131 

n.a. 

Sample:  k=0 

0.9000 

35.05 

i 

0.82 

0.9299 

13 

k=l 

0.9000 

24.92 

0.91 

0.2257 

38 

Not  appl icable 

MTwo  replications  chose  m=0  and  ninety-eight  chose  m=l  . Ihe  choice 
is  theoretically  correct. 


n.a . 
4332 
991 

n.a . 
2645 
3915 

n.a . 
718 
67b4 
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the  normal  distribution  with  mean  0 and  variance  var(o),  given 
in  {29).  Therefore  Q /var(6  j , which  is  the  basis  of  the  theoretical 
intervals  in  Table  2,  provides  centered  intervals  which  have  the  shortest 
possible  widths.  Since  the  confidence  intervals  based  on  the  sample  data 
rely  heavily  on  the  procedure  in  Section  2,  there  is  no  reason  to  expect 
them  to  be  the  shortest  possible.  Nevertheless  the  mean  width  of  .9299 
for  k=0  and  0 =.g  indicates  that,  at  least  in  this  case,  a scheme  with 
k=0  has  little  focusing  power.  We  examine  this  issue  in  greater  detail 
shortly. 

Notice  that  m = 1 for  0 =.l  and  k=l . This  agrees  with  theory  since 
every  time  a job  enters  service  immediately  upon  arrival,  a renewal  occurs. 

In  particular  the  mean  number  of  renewals  for  0 =.l  and  w=0  is  [5] 
n(l-p)  = 819.2,  which  does  not  differ  substantially  from  the  reported  991. 

In  section  3 we  discussed  a procedure  for  combining  results  for  k=0 
and  k=l  to  obtain  shorter  interval  estimates.  Table  3 lists  the  results 
based  on  (12)  and  compares  them  with  the  theoretically  shortest  achievable 
interval  for  l-a=.9,  the  smallest  achievable  theoretical  probability.  The 
dramatic  reduction  in  widths  compared  to  those  in  Table  2 is  apparent.  Notice 

Table  3 

Intersecting  Confidence  Intervals 


Coverage 

Interval 

0 

Rates 

Widths 

0.1 

Theoretical 

.90 

0.0551 

Sample 

.96 

0.0694 

O.b 

Theoretical 

.90 

0. 1918 

Sample 

.95 

0.2472 

0.9 

Theoretica 1 

.90 

0.1/88 

Sample 

.74 

0.1666 

i 

; 


H 


that  the  rates  for  o = .1  and  .5  remain.  The  poor  performance  for  o = .9  is 
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to  be  expected  since  Table  2 dictates  a maximal  achievable  rate  of  .82. 

We  next  discuss  the  poorer  than  expected  performance  for  o = .9.  Table 
4 shows  the  frequency  of  intervals  corresponding  to  the  seven  cases  enum- 
erated in  Table  1.  Notice  that  for  0 =.l  and  .5  for  k=0,l  the  intervals 
occur  principally  among  cases  1,2  and  3.  For  0 =.9  and  k=l  case  3 
occurs  exclusively.  However  for  0 =.9  and  k=0  the  less  desirable  cases 
4 and  5 occur  in  98  replications.  In  particular  the  91  of  case  5 offer 
insight  into  why  the  interval  width  for  this  case  is  so  large  in  Table  2. 


Table  4 

Empirical  Frequency  of  Interval  Estimates  by  Case  ( n=81 92 ; 


1 

5 

Interval 

ii 

o 

k=l 

k=0 

k=l 

[ava2] 

80 

98 

58 

88 

[a1 ,1] 

1/ 

0 

26 

0 

[0,a2J 

0 

1 

0 

12 

[0 , a i ]u  fa  2 , 1 J 

3 

0 

7 

0 

[0  ,a  ^ ] 

0 

0 

9 

0 

[a2,l] 

0 

1 

0 

0 

[0,1] 

0 

0 

0 

0 

.9 

k=0  I k=l 


Here  f(0)  in  (10)  is  inverted  from  the  desirable  situation  that  arises 
in  cases  1 ,2  and  3. 

Before  passing  final  judgement  it  is  instructive  to  investigate  the 
effect  of  increased  sample  size  on  coverage  rate  and  interval  width  for 
0=.9.  Table  5 compares  results  for  n=8l92  and  n=16384.  Notice  the 
substantial  improvement  in  coverage  rate  for  n=16384  and  k=0.  Regrettably 
no  similar  improvement  occurs  for  the  interval  width.  Moreover, 
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lable  5 

Comparison  of  Results  for  o = .9  and  n=8192,  16384 


Mean 

-4 

Variance  x 10 

Coverage 

Rate 

Interva 1 
Width 

Order  of 

Schane 

_★* 

m 

n=8192 

Theoretical 

0.9000 

29.54 

0.95 

0.2131 

n.a . 

Sample:  k=0 

0.9000 

35.05 

0.82 

0.9299 

13 

k=l 

0.9000 

24.92 

0.91 

0.2257 

38 

n= 16384 

Theoretica 1 

0.9000 

14.77 

0.95 

0.1507 

n.a. 

Sample:  k=0 

0.8996 

12.00 

0.95 

0.9275 

13 

k=l 

0.8996 

12.82 

0.94 

0.2028 

38 

Intersecting 

Intervals 

n=8192 

Theoretical 

0.9000 

n.a. 

>.90 

0.1788 

n.a . 

Sample 

0.9000 

n.a . 

0.74 

0.1666 

n.a . 

n= 16384 

Theoretica  1 

0.9000 

n.a . 

o 

Ai 

0.1264 

n.a . 

Sample 

0.8996 

' n.a. 

0.89 

0.1518 

n.a . 

No.  of 

Renewa 1 

N - 1 
n 


n .a . 

718 

6/64 

n.a. 

1362 

13537 


n.a. 

n.a. 

n.a. 
n.a . 


— H 


The  use  of  intersecting  intervals  offers  little  improvement  over  k=l. 

A check  of  the  interval  case  frequencies  in  Table  6 shows  that  while  case 
5 occurs  less  frequently  for  the  larger  n,  it  is  still  dominant  for 
k=0. 


Table  6 

Case  Frequency  Comparison  for  0 =.9 


Case 

Interval 

k= 

n=8192 

0 

n= 16384 

k= 

n=8192 

1 

n= 16384 

l 

[ara^] 

1 

10 

0 

U 

2 

[a  ^ * l ] 

0 

13 

0 

0 

3 

[0,a.J 

1 

0 

100 

100 

4 

[0,a,]  [a2,l] 

7 

5 

0 

0 

5 

[O.a^  ] 

91 

72 

0 

0 

6 

[a2 ,1 ] 

0 

0 

0 

0 

7 

[0,1] 

0 

0 

0 

0 

9.  Recommendations 

The  results  in  Section  8 offer  an  encouraging  picture  for  fractile 
estimation  and  provides  evidence  on  how  to  judge  computed  interval  estimates 
for  their  usefulness.  Based  on  these  results  one  can  recomnend  the  following 
steps: 

1.  Use  the  computing  schemes  in  Sections  6 and  the  criterion  (27) 

in  Section  7. 

2.  If  the  intervals  are  case  1,2  or  3 for  k=0  and  k=l,  use  these 

intervals  and,  if  desired,  form  the  intersecting  interval 
(with  lower  probability)  as  in  Section  3. 
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If  cases  4,  5,  6 or  7 arise  for  k=0  do  not  use  the  interval. 

Do  likewise  for  k=l. 

If  the  number  of  renewals  turns  out  to  be  small  do  not  use  the 

intervals  since  the  applicability  of  asymptotic  results  remains 
in  question. 
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Appendix 


The  appendix  contains  algorithms  for  computing  point  and  1-.* 
confidence  interval  estimates  for  0 = Pr{Wi  £ w)  for  specified  w,  a and 
sample  size  n. 


RUNS  converts  a sample  record  W(1 ) , . . . ,W(n)  to  a sequence  of  run  lengths 
L(1 ),. . . ,L(M)  , where  M £ n,  and  computes  a point  estimate  of  0. 

INPUT 

JMAX  h upper  bound  on  M to  be  considered. 

FRACTW  = w. 

n = number  of  observations  in  original  W record. 

{W(i )}  = sample  record. 

OUTPUT 

THETA  = point  estimate  of  0. 

JMAX  = largest  value  of  M used. 

n h actual  number  of  observations  used. 

{ L ( j ) } 5 sequence  of  run  lengths. 


KUTS  computes  the  outpoints  for  rule  m in  (27). 
INPUT 


JMAX,  n 
k = 


{ 


{L ( j ) } from  RUNS  output. 

0 if  recurrent  state  is  all  zeros, 

1 if  recurrent  state  is  all  ones. 


OUTPUT 


JM  = number  of  cutpoints  £ JMAX/2  + 1 . 

{M( j ) } 5 sequence  of  cutpoints  for  j=l,...,JM  . 

{ COUNT ( i ) > s frequency  distribution  of  run  lengths  for  k's. 

NOTE:  This  algorithm  defines  M( 0)  0 and  can  be  dropped  if  the 

zero  subscript  is  not  permitted.  The  COUNT  sequence,  a I though 
not  required,  may  assist  one  in  determining  the  degree  of  confidence 
to  have  in  the  overall  interval  estimation  procedure.  See  Section  9. 


I 
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CALC  computes  the  sample  variances  described  in  Section  6. 

INPUT 

JMAX  and  (L(j)}  from  RUNS  output.  * 
m = length  of  the  recurrent  state  y . 
k = from  KUTS  input,  k=0,l  . 

OUTPUT 

C = sum  of  the  C.  as  defined  in  (22). 

J 

S = sum  of  the  S.  as  defined  in  (22). 

J p 

SC  = sample  variance  of  the  C.,  s (C)  as  defined  in  (23). 

J o 

SS  = sample  variance  of  the  S ^ , s (S)  as  defined  in  (23). 

SCS  = sample  covariance  of  Cj  and  S^,  s(C,S)  as  defined  in  (23). 

N = N -1  as  defined  in  (22). 

NOTE:  This  algorithm  may  be  used  for  any  m between  I and  M(JM), 

where  M(JM)  may  be  found  from  the  output  of  KUTS.  Tor  m-(),  which 
is  the  assumption  of  independence,  one  should  not  use  this  routine 
but  notice  that:  C=n,  S=nTUETA,  SC=SCS=0,  N=n  and  SS  nTHETA( 1 -THETA) , 

where  THETA  is  taken  from  the  output  of  RUNS. 


Cl  calculates  the  endpoints  of  the  1 - a confidence  interval  for  o and 
estimates  the  variance  of  Z defined  in  (7). 

INPUT 

THETA  = point  estimate  for  0 from  the  output  of  RUNS. 

ZALPHA  = the  (l-a/2)  quantile  of  the  standard  normal  distribution. 
C,  S,  SC,  SS,  SCS  from  the  output  of  CALC. 

OUTPUT 

SZ  = s^(Z)  as  defined  in  (25). 

A1  = lower  limit  of  confidence  interval. 

A2  = upper  limit  of  confidence  interval. 

B2  = coefficient  of  quadratic  term  of  f(o)  in  Table  2. 

DISC  3 discriminant  for  f(O)  in  Table  2. 

NOTE:  If  DISC  is  negative  then  one  has  Case  7 of  Table  2.  If  112 

is  negative  one  has  an  inverted  interval.  Cases  4,5  or  6 of  Table  2. 


A1  gori  thin  RUNS 

I.  j < 1 . 

2-  L(j)  < 0. 

3.  j j + 1 . 

4.  If  j ± JMAX  go  to  2. 

5.  j -<-  1 . 

6.  XOLD  <-  0 . 

7.  L ( 1 ) < -1  . 

8.  i 1 . 

9.  THETA  ^ 0 . 

10.  L( j ) < L ( j )+l  . 

II . XNEW  -<  0 . 

12.  If  W(i)  <_ /RACTW  then  XNEW  < 1 . 

13.  THETA  --  THETA  + XNEW  . 

14.  j <-  j + (XNEW-XOLD)2  . 

15.  i «-  i+1  . 

16.  XOLD  «-  XNEW. 

17.  If  j 5.  JMAX  and  i n go  to  10. 

18.  JMAX  < j-1  . 

19.  n <-  i-1  . 

20.  THETA  ^ THETA/n  . 

21 . RETURN  . 


A 


(j 
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Algorithm  KUTS 

1 . i +-  0 . 

2.  i <-  i+l  . 

3.  COUNT ( i ) -f-  0 . 

4.  If  i < n go  to  2. 

5.  j <-  -1  . 

6.  j +-  j+1  . 

7.  COUNT ( L ( 2j+k+l ) ) <-  C0UNT(L(2j+k+l ) ) + 1 . 

8.  If  j < i .p. [(JMAX-k-1 )/2]  go  to  6. 

9.  M(0)  0 - 

10.  j ■*-  0 . 

11.  i -f-  0 . 

12.  i +-  i+l  . 

13.  If  COUNT ( i ) = 0 and  i < n go  to  12. 

14.  If  M(j)  = i go  to  17. 

15.  j j+l  . 

16.  M(j)  i . 

17.  j ^ j+l  - 

18.  M(j ) <-  i+l  . 

19.  If  i < n go  to  12. 

20.  JM  j-1  . 

21 . RETURN  . 


1 


f! 


i 

■ 
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Algorithm  CALC 

1 . 

J *■  . 

2. 

j «-  j+1  . 

3. 

If  L(2j+k+l ) < 

4. 

K1  j • 

5. 

j <-  i .p.  [( JMAX 

6. 

j < j-1  • 

7. 

If  L(2j+k+l ) < 

8. 

K2  < j . 

9. 

N «-  -1  • 

10. 

S «-  -m  • 

11. 

P <-  -m  • 

12. 

C <-  0 . 

13. 

LENGTH  <-  -in  . 

14. 

R <-  0 . 

15. 

Q < 0 . 

16. 

SS  <-  0. 

17. 

SC  o . 

18. 

SR  0- 

19. 

j <-  K1  - 1 . 

20. 

j «-  J+1  • 

21. 

X «-  L(2j+k+l ) . 

22. 

If  X < m go  to 

23. 

C <-  C + X + LEf 

24. 

S « S + X + P • 

25. 

R < R + Q . 

27.  SC  <-  SC  + X - m + (LENGTH  + m)2. 

28.  SS  SS  + X - m + (P  + in)2. 

29.  SR  SR  + Q2. 

30.  LENGTH  + 0 . 

31 . P «-  0 . 

32.  Q ■«-  0 . 

33.  X «-  0 . 

34.  If  j > K2  go  to  41 . 

35.  LENGTH  LENGTH  + X. 

36.  P P + X • 

37.  X «-  L(2j+k+2) . 

38.  LENGTH  •«-  LENGTH  + X. 

39.  Q Q + X . 

40.  Go  to  20. 

41 . If  N < 0 go  to  49. 

42.  SC  (SC  - C2/N)  . 

43.  SS  < (SS  - S2/N). 

44.  SR  ♦ (SR  - R2/N). 

45.  SS  4-  kSS  + (l-k)SR  . 

46.  SR  < kSR  + (l-k)SS. 

47.  S < kS  + (l-k)R. 

48.  SCS  (SC  + SS  - SR)/2. 


49.  RETURN  . 
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A1  gorithin  Cl 

1 . A1  «-  0 . 

2.  A2  ♦ 1 . 

3.  SZ  «-  SS  - 2THETA(SCS)  + THETA2SC. 

4.  B2  «-  C2  - ZALPHA2SC  . 

5.  B1  <-  2(C  S - ZALPHA2SCS). 

6.  BO  S2  - ZALPHA2SS  . 


7.  DISC  <-  B1  - 4B0  B2  . 

8.  If  B2  = 0 go  to  19. 

9.  If  DISC  < 0 go  to  18. 

10.  A1  (B1  - DISC2)/ (2B2). 

11.  A2  (Bl  + DISc'2)/(2B2)  . 

12.  If  B2  > 0 go  to  16. 

13.  TEMP  <-  Al. 

14.  Al  «-  A2. 

15.  A2  t-  TEMP. 

16.  If  Al  < 0 then  Al  «-  0 • 

17.  If  A2  > 1 then  A2  «-  1 . 

18.  RETURN  . 

19.  If  B1  > 0 then  Al  < B0/B1. 

20.  If  B1  < 0 then  A2  <■  B0/B1. 


21 .  Go  to  16. 
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For  each  m=l,2,...  the  sequence  is  demarcated  by  entrance  to  the  recurrent 
event.  The  subsequences  between  the  demarcation  points  thus  form  independent 
epochs  by  assumption.  Classical  techniques  then  yield  variance  estimates  for 
the  number  of  ones  and  zeros  in  the  epoch  as  well  as  an  estimate  of  the  co- 
variance  of  the  ones  and  zeros.  A quadratic  equation  in  o is  solved  to 
obtain  an  interval  estimate. 

Each  string  of  all  ones  or  all  zeros  examined  yields  a 1-a  confidence 
interval.  The  intervals  are  intersected  to  obtain  shorter  intervals  with 
confidence  greater  than  l-2a.  Since  each  m=l,2,...  yields  an  interval,  a 
conservative  rule  is  developed  to  determine  the  m whose  interval  is 
finally  used.  This  rule  is  based  upon  the  empirical  run  lengths  in  the 
binary  sequence. 

This  procedure  is  then  applied  to  interval  estimation  of  the  fractile  for 
the  waiting  time  distribution  in  a simulation  of  the  M/M/1  queue  with  activity 
level  0.9.  For  0=0.1  and  0.5  the  proposed  method  worked  well.  For  0=0.9 
results  showed  some  degradation.  An  error  analysis  led  to  a set  of  recom- 
mendations for  keeping  performance  in  practice  close  to  the  desired  theoretical 
levels.  An  appendix  describes  algorithms  for  computing  the  critical  quantities 
upon  which  the  proposed  method  relies. 
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